#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "matrixop.h"

#include "functions.h"

// 
void lanczos_method ( double ** A, double * x , int size, int k){

	double ** Q;	// la matrice de départ



	double * q = malloc( size * sizeof( double ));
	double * r = malloc( size * sizeof( double ));
	double * v = malloc( size * sizeof( double ));
	int i, j;
	
	// on alloue l'espace pour la matrice
	Q= malloc( size * sizeof ( double *));
		for(i = 0; i < size; ++i){
			Q[i] = malloc( size * sizeof(double));			
		}

	double alpha, beta;
	
	q = div_vect_scal( x, size, vector_norm(x, size));
	for( j = 0; j < size; ++j){
		Q[j][0] = q[j];
	}
	
	r = mul_mat_vect ( A, q, size);
	alpha = prod_vect_vect( q, r, size);
	r = sous_vect_vect( r, prod_vect_scal( q, size, alpha), size);
	beta = vector_norm(r, size);
	
	for( i = 1; i < k; ++i){
		printf(" alpha = %lf\t  beta= %lf\n",alpha, beta);
		for( j = 0; j < size; ++j){
			v[i] = q[i];
		}
		
		q = div_vect_scal( r, size, beta);
		for( j = 0; j < size; ++j){
			Q[j][i] = q[j];
		}		
		r = sous_vect_vect( mul_mat_vect ( A, q, size), prod_vect_scal( v, size, beta), size);
		
		alpha = prod_vect_vect( q, r, size);
		r = sous_vect_vect( r, prod_vect_scal( q, size, alpha), size);
		beta = vector_norm(r, size);
		
		if(beta == 0)
			printf(" beta est NULL \n ");			
	}

	printf(" ça a l'air de marcher \n");
	for(i = 0; i < size; ++i){
		
			for(j = 0; j < size; ++j){
				printf(" \nQ[%d][%d] = %lf",i,j, Q[i][j]);
				
			}
		}
}




